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Turbulent Erosion of Magnetic Flux Tubes 
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ABSTRACT 

Results from a numerical and analytical investigation of the solution of a nonlinear axially symmetric 
diffusion equation for the magnetic field are presented for the case when the nonlinear dependence of 
the diffusivity v(B) on the magnetic field satisfies basic physical requirements. We find that for 
sufficiently strong nonlinearity (i.e. for sufficiently strong reduction of v inside the tube) a current sheet 
is spontaneously formed around the tube within one diffusion timescale. This sheet propagates inwards 
with a velocity inversely proportional to the ratio of the field strength just inside the current sheet 
to the equipartition field strength Bq/B c , so the lifetime of a tube with constant internal flux density 
is increased approximately by a factor not exceeding Bq/B c , even for infinitely effective inhibition of 
turbulence inside the tube. Among the applications of these results we point out that toroidal flux 
tubes in the solar convective zone are subject to significant flux loss owing to turbulent erosion on a 
timescale of ~ 1 month, and that turbulent erosion may be responsible for the formation of a current 
sheet around a sunspot. It is further proposed that, despite the simplifying assumptions involved, our 
solutions correctly reflect the essential features of the sunspot decay process. 

Subject headings: magnetic fields — MHD — turbulence — sunspots 

1. Introduction 

The relation of turbulence to magnetic structures is ambivalent. On the one hand, any magnetic field concen- 
tration will tend to be sm eared out by the process of turbulent diffusion due to smaller scale motions. Since first 
proposed by Sweet (1950), this process of turbulent magnetic diffusion characterized by a diffusivity v ~ Iv (where 
I is the scale of the turbulence and v is its r.m.s. velocity) was and, despite repeated controversies, has remained 
a cornerstone of mean field magnetohydrodynamics, observable e.g. in the form of a slow dissolution of unipolar 
areas on the solar surface ( Sheeley, DeVore, fc Boris 1985 ). On the other hand, turbulent motions tend to sweep 
the field lines together at convergence points of the flow component perpendicular to them; thus, the field becomes 
concentrated in flux tubes (Weiss, Proctor fc Weiss 1964, 1982). This process of flux concentration is thought to 



be aided by the fact that magnetic fields tend to inhibit turbulent motions, and thus turbulent diffusion. 

In homogeneous stationary MHD turbulence, the fine structure of the magnetic field is essentially determined 
by a dynamical equilibrium of these two main processes. This leads to an intermittent field structure where at 
any given instant of time most of the magnetic flux is concentrated into thin flux tubes of about equipartition 
strength, which however possess a lifetime not exceeding by much the turbulent correlation time. In inhomogeneous 
and/or non-stationary situations, however, it may often happen that a flux tube gets into an environment where 
its flux and/or field strength far exceed the values corresponding to the local turbulent equilibrium described 
above. In such situations one may expect that turbulent diffusion will dominate over concentration, leading to 
the slow dissolution of the tube. A case in point are the sunspots, which are good examples of such oversized 
and over-intense flux tubes (as compared to the local turbulence scales). They are seen to decay in a few weeks 
after their formation in a way very much reminiscent of turbulent diffusion (Zwaan 1968). A similar process of 
dissolution might be expected in the case of the toroidal flux bundles lying near the bottom of the solar convective 
zone: they may provide another example of oversized and over-intense flux tubes if their sizes and strengths are 
inferred correctly from flux loop emergence models. A detailed analysis of such a turbulent decay of magnetic flux 
tubes is however complicated by the fact that magnetic fields with energy densities comparable to or exceeding 
the turbulent kinetic energy density inhibit turbulence, and thereby reduce v by a significant amount. The exact 
form of the relation v{B) is poorly known. (See Kitchatinov, Pipin, & Riidiger 1994 for a formula valid under 
certain assumptions.) However, from energy considerations it is clear that in the limit — > it should tend to 
the kinematic value uq ~ Iv, while, in the strong field limit B — > oo, it must go to zero. At the equipartition field 
value one may expect a reduction by a factor of order unity. 

The strong inhibition of turbulence inside tubes with field strengths well above the equipartition value (as in 
the case of the two examples mentioned above) suggests that the decay mostly advances by the "gnawing" action 
of the external turbulent motions on the edges of the tube. This was first proposed by Simon & Lcighton (1964), 



who called this process the erosion of the flux tube. They proposed that erosion by supergranular motions was 
responsible for the decay of sunspots. Smaller-scale granular motions may however be better suited for this task 
than the large-scale supergranules. 

This paper is an attempt at the detailed study of the turbulent erosion process of magnetic flux tubes out- 
lined above. For simplicity we consider an axially symmetric magnetic flux tube and disregard the variation of 
physical quantities along the axis of the tube, assuming that their lcngthscale is large compared with the tube 
radius. For the nonlinear dependence of the magnetic diffusivity v on the field strength B we consider a simple 
analytic function satisfying the basic physical requirements. We note that the magnetic field in general leads 
to an anisotropic diffusivity, so that, strictly speaking, our v is just the rr-component of the diffusivity tensor 
in cylindrical coordinates. With these assumptions we write down and numerically solve the nonlinear diffusion 
equation governing the problem in § |^. We show that a steep diffusion front is easily formed at the boundary of 
the tube when the internal field strength is large compared to the equipartition value. This front slowly advances 
into the tube, thus removing magnetic flux from it that is diffused outward into the surrounding weak field region, 
contains an analytical study of the problem, including the calculation of form-invariant, propagating solutions 



(§ 3.l|) and estimates for the velocity of the front (§ 3.2) and of its thickness (§ 3.3). In § W we will consider some 



implications of our results for toroidal flux tubes rising in the convection zone and sunspots. 

The problem of diffusion with a diffusivity dependent on the variable being transported appears also in other 



physical and astrophysical contexts such as turbulence propagation or viscous gravity currents (Barenblatt 1983 



Gratton & Minotti 199C). An interesting case is heat conduction in a plasma with conductivity, k, given by 



Spitzcr's formula (ft cx T 5 / 2 ). The strong dependence on T gives rise in this case to conduction fronts, which 
are relatively abrupt transitions from high to low temperatures. In that case, the front becomes steeper toward 
low temperatures. In the magnetic case, in contrast, the diffusivity is suppressed by the higher fields. Hence, we 
expect the magnetic diffusion front (i.e. the current sheet) to become steeper toward higher field intensities, as 
in a mesa-like mountain profile (cf. Fig. Id below). An important difference between the thermal and magnetic 
problems is the possibility of a stationary situation: in the thermal problem one can compensate for the energy 
being conducted with adequate amounts of heating and cooling on the hot and cold ends, so that a stationary 
situation may ensue (as found in different astrophysical contexts). In the problem of a magnetic diffusion front 
around a flux tube, in contrast, there are no source terms for the magnetic flux inside the tube, so stationarity 
can only be achieved, if at all, through an external advection term which brings back to the tube the magnetic 
flux that the ohmic diffusion is trying to remove from it. 



2. The nonlinear diffusion problem 

2.1. Formulation of the problem 

As a first approach to the problem of the nonlinear diffusion of magnetic field from a magnetic flux tube, we 
set out to solve a nonlinear one-dimensional diffusion equation with axial symmetry. The one-dimensionality may 
limit its application to actual solar problems (see Sect. but it has the advantage of retaining the main nonlinear 
effects while keeping the problem relatively simple. The diffusion equation for the magnetic field in cylindrical 
coordinates reads 



CD 



where the notations are r for radius, t for time, B for magnetic flux density and v for magnetic diffusivity. The 
inhibiting action of the magnetic field on turbulence is contained in the form of the function v(B) which should 
therefore have the following properties: 

(a) lim B ^ v = vq where vq is the unperturbed value of the diffusivity 

(b) v/vq £ 0(1) and (1 — v/vq) E 0(1) if B ~ B , where B is that value of the field intensity for which v is 

reduced by 50%. We can expect B e to be order of the equipartition flux density corresponding to a magnetic 
energy density that equals the kinetic energy density of turbulence. 

(c) v/vq -C 1 for B/B c ^> 1 (assuming that the molecular diffusivity is negligible) 
A simple function satisfying the above requirements is 

= mw: (2) 



The parameter a v determines the steepness of the diffusivity cutoff near B c . All the results presented below were 
computed with the above form of v(B); however, test runs with other forms of v(B) compatible with conditions 



(a) through (c) yield qualitatively similar findings. Furthermore, in § 3.2 below we will argue that the general 
behavior of the solutions should be similar for all forms of v(B) obeying the criteria listed above. 
As initial condition we again choose simple analytic B(r) profiles of the form 



B(r) 



Bo 



1 + (r/ro)"- 8 

with as = 22 (to model a tube with constant field strength inside), or 

B(r) = B Q exp(~r 2 /r 2 ) 
for a tube with a field profile gradually decreasing outwards. 



(3) 



(4) 



2.2. Numerical solutions 

Equation (|l|) is a nonlinear flux-conserving equation that may be solved by standard finite-differencing tech- 
niques. A two-step Lax-Wendroff scheme was applied for this purpose. In the numerical calculations, and whenever 
dimensionless quantities are implied in the paper, the following units were used: ro = Vq = B e = 1. (Note that 
the unit of time is thus the diffusive timescale Tq/vq.) 

For the solution one should also specify the boundary conditions. At r = the boundary condition is set by 
the requirement of symmetry, while at the outer boundary we experimented with different forms of the boundary 
condition to find that if this boundary is sufficiently far away (at r — 10) the form of the boundary condition exerts 
a negligible influence on the results in the physically interesting regime r < 2.5 (in the figures, only this interior 
regime is shown). Consequently, for all calculations we set the outer boundary at r = 10 as dB/dr = 0. The 
outer fictive points necessary for the Lax-Wendroff method were computed with the neglect of third derivatives; 
this method preserves the second-order accuracy of the scheme. 





Fig. 1. — Nonlinear turbulent diffusion of a magnetic flux tube. Shown are snapshots of the magnetic field profiles at 
10 equidistant instants of time (with timesteps Ato, in units of the diffusive timescale) for a magnetic diffusivity of the form 
(2) with ol v — 0, Ato = 0.5 (linear limit, panel a), a v — 2, Ato = 1 (b), a v = 3, Ato = 1 (c), and ol v — 7, Ato = 0.25 (d). 
Nondimensional units as defined in Section 2.2 

The stability criterion At < (Ax) 2 /2u seriously limits the timestep; therefore, owing to the finite available 
CPU time, the maximal spatial resolution that could be attained was 8000 grid points (uniformly distributed in 
the range < r < 10), with a corresponding timestep of 10~ 5 . 



Figure |l] presents the time evolution of the field profiles with initial condition given by (|3|) with a B — 22 and 
for different values of the coefficient a v in equation (||), namely a u — (fully linear case, panel a) and other 
three which yield increasingly steep profiles (panels b, c and d, with a„ = 2,3,7, respectively). For low a v the 
solution is still similar to the linear constant-diffusivity case, viz., it is characterized by a monotonic decrease 
of field gradients, smoothing out the field distribution, and by a gradual decrease of the flux density at r = 0. 
However, for higher values of a v a qualitatively new behaviour sets in: in a short interval of radius the field 
gradient (or equivalently the current density = —Ait/c dB/dr) greatly increases, i.e. a current sheet is formed 
spontaneously. The flux density on the axis remains constant until the arrival of the current sheet and then it 
suddenly drops to nearly zero. The sheet moves inwards with an essentially constant speed, thereby leading to 
a nearly parabolic decay law for the flux and area inside it (cf. Fig. ||). Within about the equipartition radius, 
r < r c , where B(r c ) = B e , the overall form of B(r,t) is consistent with a dependence of the form B(u) where u 
is the distance from the current sheet, with no explicit time dependence. In what follows we will call solutions 
of this type erosional solutions, a name which reflects the fact that the formation and inward propagation of the 
current sheet is due to the eroding action of turbulent magnetic diffusivity on the outskirts of the tube, while in 
the tube interior diffusion is practically absent. 
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Fig. 2. — The area enclosed within the current sheet as a function of time. The cases a v — 3 (solid) and a v = 7 
(dashed) are shown (corresponding to panels c and d in Fig. 1). 

The thickness of the current sheet quickly decreases with increasing values of a v . This limits the possibility 
of calculating numerically the solutions for high a v . The highest value for which we still had good resolution, 
a v = 7, was reached by keeping Bq at a moderate value of 3 (Fig. |T](f). In this case, the current sheet extends just 
to B ~ B c ; from there outward there is a range of radii where the solution quickly approaches the solution of a 
linear diffusion problem. For higher values of a u , this transition range shrinks to basically just the point where 
B = B c . In that case a near discontinuity in the derivative of the field profile (i.e., in the current density) appears. 

Figures |^ and |] show the propagation velocity w of the current sheet as a function of a v and B , respectively. 
The velocity w was determined by measuring the horizontal separation of consecutive B(r) curves at different 
times and B values inside the current sheet. For a constant internal field strength, l/w is a good estimate of the 
total lifetime of the tube. While the lifetime increases approximately linearly with Bq, it seems to saturate to a 
finite asymptotic value with increasing a v when Bq is constant. This asymptotic value is approximately Bq. 

The magnetic flux that is being removed from the strong-field region is diffused away outward from the current 
sheet. The rate of removal of flux is roughly 2irrcswB , where res is the radius of the current sheet. Given 
the form invariance of the solutions shown in the panels of Figure [I] with higher a„'s, we can expect the flow of 
magnetic flux across the current sheet to be essentially constant: flux cannot accumulate (nor the contrary) within 
the sheet if it is to keep a constant shape in time. Hence, the following quantity, the magnetic flux removal rate, 
[wB — v{B)dB / dr] should be constant in the whole strong-field range and across the current sheet. This is borne 
out by the numerical solutions: Figure || shows that the profile of that quantity (long-dashed curve) is basically 
constant from r = out to the neighborhood of B = B c . Thereafter it smoothly increases: in the diffusive part 
there cannot be a form invariance, since the magnetic flux must increase along time to accommodate what is being 
removed from inside. 

The existence of this class of erosional solutions is not limited to the case of the initial conditions given by 
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Fig. 3. — Current sheet velocity w as a function of a v for models with Bo = 10, as = 22 with the limits given by eqs. 
(16) [solid] and (15) without the 3~ 1,/2 factor [dashed] 




Fig. 4. — Inverse of current sheet velocity 1/ro as a function of Bo for models with a„ = 3, as = 22, with the limits 
given by eqs. (16) [solid] and (15) without the 3 -1 ^ 2 factor [dashed]. 
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Fig. 5. — Plot of w(B — B ) — v(B)dB / dr (long-dashed curve) for an intermediate time of the solutions of Fig. 1 
with ct v = 3. In the figure, the numerical solution for B (solid line) and the corresponding analytical propagating solution 
(short dashes) are shown with ordinates on the left y-axis. Two further curves are shown with ordinates on the right y-axis: 
—u{B)dB/dr (dash-dot) and —w(B — B ) (dash-triple-dot). The constancy of the long-dashed curve to the left of and 
within the current sheet is a direct consequence of the form-invariance of the erosional solutions. 



equation (|3|) . Figure |6| shows that such solutions are also found with sufficiently high values of a u for other initial 
conditions, the only difference being that for an internal magnetic field depending on r the propagation velocity 
of the current sheet is not constant anymore. 




Fig. 6. — Same as in Fig. 1 but with a Gaussian profile as initial condition [eq. (4)], oi v = 3.5, Ato = 1. The velocity 
of erosion is no longer constant, but, rather, decreases for increasing field on the strong end of the current sheet. 

Up to this point we tacitly assumed that the value uq of the diffusivity outside the tube is constant. This 
corresponds to a case where the size of turbulent eddies is small compared to the size of the tube. In the 
alternative case when the turbulent eddies are large, only the small-scale turbulence (with scales below the tube 
radius) may contribute to the diffusivity, while larger-scale motions will simply lead to the translational motion 
of the tube as a whole. Assuming a Kolmogorovian spectrum, the diffusivity scales as rcs 4 ^ 3 - ^ solution with 
such a rescaling of Pq in each timestep to the appropriate value of res (assuming vq = 1 at t = 0) is illustrated 
in Figure ^. The decay of the tube suffers a strong deceleration with time, making the total lifetime of the tube 
effectively infinite; however, a significant fraction (e.g. 50 or 90 percent) of the total flux is still lost on a timescale 
comparable to that of the solution without rescaling. 
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Fig. 7. — Same as in Fig. lc but using a cutoff for the size of the turbulent elements equal to the tube radius and a 
Kolmogorov spectrum for the turbulence, with Ato = 2. The value of vo in equation (2) then becomes a function of the 
tube radius and the decay of the tube is strongly decelerated 



3. Analytical solutions and estimates 

Among the features of the numerical solutions discussed above, some merit particular attention. In the present 
subsection we clarify and substantiate some of them using analytical considerations. This allows us as well to 
extrapolate them to the numerically unattainable limit a u — > oo. 



3.1. Propagating solutions 

The numerical results show that in the case of a constant internal flux density, for sufficiently strong nonlincarity 
(i.e. for sufficiently high values of a, v and/or Bq), the solution in the regime B ^> B c approximates a form-invariant 
propagating solution of type B = B(r + wt), with w approximately constant. This invariant shape is reached 
within a diffusion timescale Tq/vq. 

This prompts the question of the existence of purely propagating, form-invariant solutions of the nonlinear 
diffusion equation ([!]) . To obtain those, we go to the cartesian limit of the equation (with spatial coordinate x) 

dcf 

and try to obtain solutions which depend on t and x through the combination u = x + wt, with w a constant 
value (a positive w means a solution moving toward smaller spatial coordinate, as in the numerical solutions). The 
corresponding equation is 

, ,dB~ 



du 



-wB 



du 



= 



(5) 



This equation has as immediate solution the following class of functions: 



v{B) 



K + wB 



dB , 



(6) 



with u c and K two integration constants. Equation (^]) gives B implicitly as a function of u — u c . This class 
of solutions (II) has a flat, exponential asymptote for B — > —K/w, which, therefore, may serve to describe the 



flat region shown by the numerical solutions of the previous section toward r — > 0. We then write —K/w 
introduce the notation b = B / B e , b = Bq/B b , and use the form (§) for v so that the solution (||) becomes 

db 



wj (i + i&ho(i - &/& ) 



B n 



(7) 



The integral in equation (Q) can be easily carried out analytically for integer ol v (although the solution becomes 
somewhat unmanageable for a u > 5) or, numerically, for arbitrary a v . For instance, for a v = 2 one obtains 



w 1 + b n 



b atan& - 



log 



(1 + & 2 ) 1/2 



1 - b/b 



(8) 



The form of the solutions (0) is illustrated in Figure || for four different cases, namely a v — 2 (dotted), a u = 3 
(dot-dashed), ol v — 7 (dashed) and a very steep case, a v = 50 (solid), which illustrates well the limit a v — > oo. 
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Fig. 8. — Analytic solutions for a plane parallel form-invariant propagating solution of the nonlinear diffusion equation 
(6) with Bo = 3 and a„ = 2 (dotted), a v = 3 (dot-dashed), ol v — 7 (dashed) and a v = 50 (solid), u is dimensionless in this 
figure, unit is fo/w. The near discontinuity of the current density at the lower end of the current sheet is clearly discernible 
in the steepest case. 



The increasing nonlinearity makes the transition between the diffusive and non-diffusive regions increasingly steep, 
similarly to the case of heat conduction with nonlinear conductivity coefficient. For a v = 50, in particular, the 
current sheet is very narrow, with a sharp transition to the asymptotic horizontal regime. For B < £? c , the four 
solutions are qualitatively similar and approach a linear diffusive solution. In these plane parallel solutions, there 
is no natural lengthscale imposed by the equations: solutions exist for arbitrary w > values, independently of 
Bq. A change in w does not change the shape of the solution, but only its lengthscale, which is vq/w. In contrast, 
in the actual axisymmetric proble m a length scale is set by the radius tq of the initial tube and w can no longer 
be specified arbitrarily (cf. Section 3.2 below). 




Fig. 9. — Comparison of the analytical (dashed) and numerical (solid) solutions for Bo = 3 and ol v = 7. The analytical 
solution was computed with the w value determined from the numerical solution, and the arbitrary horizontal offset was 
set to optimize the agreement. 

We expect the propagating solutions to fit the numerical solutions of the previous section increasingly well the 
higher the value of a v and b and, at any rate, for b > 1. The goodness of the fit can be seen in Figure [| which 
contains the numerical (dotted) and analytical (dashed) solutions for the case a v — 7 with Bq = 3. The fit is 
excellent in the whole non-linear region. 

The propagating solutions cannot approximate the numerical solution well outside the current sheet for fun- 



damental reasons: the current sheet is a region where magnetic flux is being taken away from the tube and 
transported outwards. To have a full form-invariant solution, there should be a mechanism to exactly dispose of 
this flux in the outer regions of the domain. The latter can be, for example, a flux-shedding boundary condition 
or some internal device that causes the loss of magnetic flux (which should play a role like the cooling in a thermal 
problem). We have none of these, and, as a consequence: (a) the numerical solution is not form-invariant in the 
weak-field range and (b) the analytical solution goes to negative field values at some distance of the flux tube (see 
Fig. 8). 

The establishment of a form-invariant propagating region along the time evolution of the system may then be 
qualitatively understood as follows. The high diffusivity in the outskirts (B <; B c ) of the tube removes the flux in 
that domain in about one diffusive timescale r^/v®, while the strong internal field at r < res effectively inhibits 
any significant flux loss from the internal region over such short timescales. As a consequence, a steep field gradient 
(i.e. a current sheet) forms. The gradient increases until the smoothing effect resulting from the diffusive flux loss 
can balance the non-linear steepening. This is exactly the balance responsible for the form-invariant nature of 
the exact analytic solution of equation (^) , which is therefore a good approximation of the solution in the regime 
B > B c . This property of the erosional solution may be utilized in deriving a formula for the thickness of the 
current sheet (cf. Section |3~3| below). 



3.2. Current sheet velocity 

From the numerical solutions, the velocity w of the inwards propagation of the solution (and of the current 
sheet in particular) appears to increase linearly with Bq and to become asymptotically independent of a v for high 
a„ values, i.e. for increasingly effective magnetic inhibition of turbulence. This implies that the lifetime of flux 
tubes remains finite even in the limit of a v — > oo, i.e. for infinitely effective inhibition of turbulence inside the 
tube! Although Fig. || strongly suggests such a behaviour, a strict proof of this result clearly cannot be given on 
the basis of numerical results alone. Below we demonstrate with approximate analytical methods that a lower 
limit for w indeed exists. Here we concentrate on the case of a flux tube with constant internal flux density and a 
corresponding form- invariant propagating B(u) profile in the erosional solution. It is however clear that for tubes 
with the same value of Bq but outwards decreasing field strength the efficiency of the erosion process can only 
increase. 



The analytical results for the plane parallel case (cf. Section 3.1 above) are of no help in the present problem 
since these solutions exist for any value of w. It is clear from a dimensional analysis that the most generic combi- 
nation of the parameters entering the problem (z^o, a„, ro, B ei Bq) with velocity dimensions is f(Bo/B c , a^vo/ro- 
Our task consists of determining the form of the function f(Bo/B c , a v ) in the limit a v — > oo. 

To obtain that function, note that in that limit the current sheet becomes infinitely narrow, is placed precisely 
at the equipartition radius, r = r c , and coincides in its full extent with the current sheet of the corresponding 
form-invariant propagating solution. Moreover, the transition between the current sheet and the v = vq regime 
occurs sharply at B = B c . In fact, v{B) will have a near discontinuity at B — B c and so will the derivative dB/dr. 
We saw in the previous section that the form invariance implied that the flow of magnetic flux, wB — v{B)^, 
was constant in the whole form-invariant section of the solution down to the lower end of the current sheet (and 
slowly departs from the constant value when going outward from it). In our case, this means that that quantity is 
constant down to B = B c . Using all these facts, we can now obtain a first estimate for the velocity of the current 
sheet. Equating the values of that quantity on both ends of the current sheet, we obtain 



' dB s 

wB = wB c -v [ — ) (9) 

ar I B=B C 



where the last derivative is meant as the limit value when one approaches B = B e from outside the current sheet. 
This leads to 

B c 

w =—~5 ( 10 ) 

n B - B e 



with 

1 def 1 fdB 



n B e V dr , B=Bc 



(11) 



The question is now how to determine r% in equation (|l0|). An upper limit can be set by noting that outside the 
current sheet, i.e. for r > r c , d 2 B(r)/dr 2 > 0, so that the straight line B = B c [1 — (r — r )/ri], tangential to the 
solution at r c , certainly runs below the actual field distribution. The latter, then, contains more magnetic flux 
than the straight line. Using the approximation <£>o — ^qBq where $0 is the total magnetic flux in the whole 



domain, we have 

/•re+ri 

2ttS c / [l-(r-r e )/ ri ]rdr < ttB^ - r c 2 ) . (12) 

J r c 

This inequality should be valid at any stage after the moment when the form-invariant propagating solution finally 
sets in and ri settles to a basically constant value. In order to determine the constraint on rj it suffices to know 
the value of r c at that instant. A lower limit for this value is obviously zero, so solving the inequality ([l2]) for r% 
we get 



3 

ri<--r e 



^ 2 i o I 2 2\ 

s r. +3-(r -r e ) 



1/2 /Q R X 1/2 

y 2 , o ^0 2 



Now, in the high-a,, limit, the current sheet is very thin, so that r e < tq. This yields 



r V4 B, 

Using ( |l0| ) we obtain the following lower bound for w in the limit B ^> B c 

-, / D \ 3/2 
2/ 1 / B c X 



ri <f? + 3^V /2 . (14) 



™>— -7= -5M ( 15 ) 

We note without going into similar details that replacing the above straight line with a logarithmic profile 
[stationary solution of the linear diffusivity equation with the boundary condition B(r e ) = B e ] would lead to a 
transcendental equation instead of the second-order equation (|l2|). Its solution yields a condition nearly identical 
to ([l5|) without the factor How good this lower bound for w is can be seen in Figures and ||, where the 
dashed line represents the right-hand-side of equation ( |l5| ) (or its inverse, in the case of Fig. ||), without the V3 
factor. 

A more realistic estimate of the value of r e at the time when the erosional solution sets in can be given taking 
into account that the initial transient phase physically consists of the period when the diffusive gnawing of the 
outer parts of the tube increases the gradient further in. Thus, one may expect the duration of the transient phase 
to be i t r ~ r o 2 /^0: i- e - the diffusion timescale. With this, one has approximately r c (t — t tr ) ~ r (l — wr /v ). 
Inserting this condition into the inequality ( fl5| ) (or, better still, into the equivalent inequality for the logarithmic 
profile discussed in the previous paragraph), one obtains a relation between w and r\. This relation, coupled with 
equation (|l0|), can be solved to find that in the limit Bq ^> B c 

^ 2 -l/3|e^ (16) 
B r 

This new lower bound (represented as a solid line in Figs. |^ and ^) turns out to give a very good approximation to 
the numerical results for w. In fact, the value of the coefficient in equation (|l^) shows only weak sensitivity 

to arbitrary order-of-unity factors in the expression of t tr and in the upper limit used in the integral constraint 
(typical changes are by ±0.1). Equations ([l6]) and ( |ic| ) then clearly also imply that r\ and ro are of the same 
order of magnitude. 

Closing this subsection we should stress that the validity of the results derived here is not confined to the 
particular choice of the v(B) function given in equation (^), as this particular form was nowhere used. For more 
general profiles, the parameter a u may simply be defined so that the relation 



vq \dB 

remains valid, i.e. a u simply parameterizes the steepness of the diffusivity cutoff at B c . It is easy to see that all 
the above results remain valid in this more general case. 



3.3. Current sheet thickness 



Both in the numerical results of Sect. 2.2 and in the propagating solutions of Sect. p.l| the thickness of the current 
sheet rapidly decreases for increasing a v . In this section we obtain an analytical estimate for this dependence and 
show that, in fact, the current sheet thickness has a power-law dependence on the internal field strength, B /B c , 
with exponent — a v . We show this for the analytical solutions of Sect. 3.1. The result is then immediately 



applicable to the numerical erosional solutions with high a u , given the excellent agreement between them. 



In the current sheet, the field distribution must go through an inflection point. Define the thickness of the 
current sheet as the scalelength of B(r) at that point. Using equation (||), we get 

where the index i refers to values at the inflection point. Substituting here Jl6|), 

* Stfia.fH)' *(*)-. ( 19 , 
ro \voj B e \B e J 

Next we show that in the high-a^ limit the inflection point occurs high up in the current sheet, i.e., B\ — ► Bq. 
At the inflection point one has — 0; using the general integral [eq. (0)] of the plane-parallel form-invariant 
propagating solution, this implies 

VidB = Bi-B ' ^ 
Substituting here the expression (^) of v{B) we find 



u \Bi J \B 

the high-a„ limit of which is indeed 



A - - B Q . (22) 



Comparing this with equation (|19|), and using the limit -B >> B c , we finally obtain: 
which was to be shown. 



-S^aJ^-) 1 ^ ■ (23) 



4. Applications 

The solar magnetic field appears in the form of magnetic flux tubes of different sizes in the solar atmosphere. 
Through the studies of the emergence of the active region fields as well as from general theoretical results of 
magnetoconvection and MHD turbulence, one assumes that also the magnetic field in the convection zone may be 
in the form of individual flux tubes or bunches of them. For the application of the one-dimensional results of the 
foregoing sections to the Sun, the magnetic flux tubes must be, even if only approximately, cylindrically symmetric. 
In our case, a necessary condition for this is that the lengthscale of variation of properties along the tube axis 
be much larger than transversely to it (e.g., the tube radius must be much smaller than all local stratification 
scaleheights) . Under this condition [typical, for instance, of the so-called thin-flux-tube (tft) approximation], 
equation ([[]) can be used to describe approximately the diffusion of field through the surrounding turbulence. 
We stress that while the assumption of cylindrical symmetry is a common feature of our models and of the tft 
formalism, yet here we break away from the tft approximation by explicitly considering the field variation across 
the tube. In the following we consider two cases of current interest, first one in which the cylindrical shape of the 
tube is a good approximation, namely, the magnetic tubes in the deep convection zone, and, second, the diffusion 
of the field in sunspots. Sunspots can be simplified using axially symmetric models; yet, they do not possess 
cylindrical symmetry, so this case can only be seen as a simple approximation to a far more complex problem. 



4.1. Toroidal flux tubes at the bottom of the solar convective zone 



It is generally accepted that most of the magnetic flux responsible for solar activity phenomena in gene ral, and 
for the appear ance of active regions in particular, is stored close to the bottom of the solar convective zone ( Bpicgc] 
fc Weiss |l98C ; see also Morcno-Inscrtis, Schiisslcr, fc Fcrriz-Mas 1992| ) . The question arises whether these toroidal 
flux bundles are affected by the phenomenon of turbulent erosion modelled above to any significant degree. Using 
equation (Q6h, the flux loss rate from a tube is 



- $ = 2irr cs B w = 2 2 / 3 7 Wo B c r cs /r - 2 2 / : W 5 o ($/$ ) 1 / 2 



(24) 



where we used the approximation res — ( < J ) / 7r -E'o) 1 ^ 2 : and we took into account that res and 4> are functions of 
time, while w is determined by the initial radius ro (or equivalently $o) i n the case when the tube is not moving 
compared to its surroundings. In a situation where a tube is moving through the surrounding plasma (as in the 
case of emerging flux loops) with a significant speed (V > vo/tq), the external flow relative to the tube will sweep 
away the weak outer field thereby continuously "reinitializing" the decay; in that case, the actual, time-dependent 
value of res should be used in the expression for w, and consequently equation (0) must be replaced by: 



- $ = 2 2 I z -kvqB c (moving tube) 



for magnetic tubes in the deep convection zone. 

2 , 



(25) 
The full 



Consider first the implications of equation 
turbulent diffusivity there is of order i/qq ~ WH <~ lOOOkm^/s, (W ~ 60m/s is the r.m.s. vertical velocity of 
turbulence), but as the scale of the turbulent eddies is the pressure scale height H ~ 5 • 10 4 km, i.e. much larger 
than the tube size, the diffusivity should be rescaled as 
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Writing this into equation (pij) we find 



$ = -WB C 
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(26) 



(27) 



Numerical simulations of the emergence of buoyant magnetic flux loops through the solar convective zone 



1986J, |1993| , |1995| , |1994| show that in order to get a satisfactory agreement with the observations of solar active 
regions, the toroidal flux tubes on which the perturbations develop must have magnetic flux densities of Bq ~ 10 5 G, 
while the equipartition field at the bottom of the convective zone is B e ~ 10 4 G. Those models do not take into 
account any diffusivity, so, to match the observed fluxes in active regions, the tubes must have magnetic fluxes of 
order 10 22 Mx also in the deep convection zone. With these values, using equation (|27| ) we find that — $/$ ~ 1 
month for the flux loss timescale. This is very much shorter than the length of the solar activity cycle. While the 
timescale increases with time owing to the factor (<&/ <&o) 1//2 m equation (|27|), the greater part of the original flux 
will still be lost on essentially the initial timescale. 

On the basis of such short lifetimes, it appears that flux tubes with the above properties cannot be stored 
in the lower part of the convective zone for times comparable to the solar cycle. We note however that and 
consequently the timescale, depend rather sensitively (as W 2 , since B e cx W) on the r.m.s. turbulent velocity. 
Thus, reducing W by a factor of 5-10 could increase the decay timescale to several years, comparable to the solar 
cycle. 

However, a careful analysis of the possible sources of error, taking into account the results of the relevant 



numerical simulations (Chan & Sofia 1989, Lydon, Fox, & Sofia 1992, Kim et al. 1996) shows that the W values 



derived here are subject to uncertainties not exceeding ±30-40 %. In particular, there seems to be no obvious way 
to reduce W by more than about 20 % at the base of the convective zone. So the flux storage problems mentioned 
above appear to be unavoidable. 

Thus, the calibration of local convection theories to numerical experiments leaves now little room to change 
W in the unstably stratified part of the convective zone. On the other hand, our present ignorance regarding the 
detailed structure of the overshooting layer below still allows the possibility of storing the toroidal flux there. 

Turning now to expression (p5|), we consider the loss of flux in unstable magnetic tubes that rise across the 
convection zone to produce active regions. To that end, we substitute for v using equation (|2l) to get: 



-W B r 



4tt$ 2 



1/3 



(28) 



A first impression for the orders of magnitude predicted by this formula can be obtained by substituting in it 
a few time profiles z(t) and B (t) typical of the tft numerical simulations of the rise of magnetic tubes mentioned 
above, with z the depth below the photosphere. Taking, for instance, for z(t) and B a (t) the position and field 
strength of the summit of the rising loop in the simulations of Caligari, Morcno-Inscrtis, & Schiissler (1995), one 
can obtain $ from equation fl28|). To avoid the added complication associated with the sliding motion of the tube's 
mass elements along the field lines, we use here time profiles corresponding to a non-rotating case (tt = 0), for 
which the mass element at the tube summit is fixed along time. On the other hand, one has to substitute for W, H 
and B c using a stratification model and the following formula?: 



(29) 



1996 , where V a d was simply taken to be 0.4 (a good approximation except near the surface) and 

Bo = AirpxW 2 , 

with the numerical factor x = 1 + 2 • 0.61 2 taken from Chan fc Sofia ( 1989| ). 



(30) 
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Fig. 10. — Flux loss from flux loops rising through the solar convective zone. Solid: $ = 10 22 Mx, Boo = 1-2 ■ 10 5 G. 
Dashed: $ = 10 23 Mx, B o = 1.2 ■ 10 5 G. Dash-dotted: $ = 10 22 Mx, B 00 = 2 ■ 10 5 G. 

Figure [l^ shows the decrease of the magnetic flux during the passage through the convective zone predicted by 
Eq. (^) for tubes of different sizes and field strengths. We first use the z(t) profiles obtained from the thin-flux- 
tube numerical simulation of the rise of tubes with 10 22 Mx and Boo = 1-2 • 10 5 G (solid curve) or Bqq = 2 • 10 5 
G (dash-dotted). We integrate Eq. ( pq ) downward starting from the surface with that value of flux. We see that 
these tubes must lose a substantial fraction of the their original magnetic flux along the rise (in contrast to the 
simplifying assumption systematically made in the thin flux tube simulations). The difference between the two 
curves is due to the shorter timescale of the tube with higher initial field strength. Anyway, both curves show flux 
values of order 10 23 Mx in deep levels. So, we repeat the calculation, this time taking z(t) from a tft simulation 
with 10 23 Mx and Bqo = 1.2 • 10 5 G (dashed curve) (but start the integration of Eq. ( pS| ) at the surface with the 
same $ value as for the other curves, for easiness of comparison). Here too the tube suffers a substantial magnetic 
flux loss. However, the calculation shows only a relatively small sensitivity to the initial value of the magnetic 
flux in the emergence computation. This is because the magnetic flux of the tube only enters the tft calculation 
via the drag force opposed by the surrounding medium to the advance of the tube. In the range of fluxes of the 
figure, the radius of the cross section is large enough for the drag force not to play a primary role in determining 
the speed of rise - hence the time evolution profile of all quantities in the tube is fairly independent of its magnetic 
flux. This is a fortunate circumstance as it justifies a posteriori our (strictly speaking inconsistent) approach of 
studying flux loss using models computed with the assumption of a constant flux. It is perhaps also adequate to 
mention at this point that the use of z(t) and B (t) profiles taken from numerical simulations for a rotating case 
(fi = Qq) yields results basically coinciding with those shown in Figure [l(] within a factor 2 at most (yet, the 
interpretation of the results for the rotating case has to be done with some care, as indicated above). 

From all this we conclude that the toroidal flux tubes whose emergence leads to the formation of solar active 
regions must lose a significant part of their flux by turbulent erosion during their passage through the convective 
zone. In fact, it would be interesting to include the formulas developed in this section into a tft model, so that one 
could calculate the time evolution of the magnetic flux in a more self-consistent way. Doing that could result in 
a non-constancy of the magnetic flux of the tube as one moves along its axis at any given instant of time. This 
would be no violation of the solenoidality of the B field, but, rather, just reflect the fact that the external sheets 
of a stretch of the tube may at some level detach from its concentrated core and go into a weak-field phase that 
is no longer counted as part of the tube in the tft model. 



4.2. Sunspots 

The sharpness of the outer penumbral boundary of sunspots (i.e. the sudden change from filamentary penumbral 
fine structure to nearly normal granulation in less than a resolution element, Zwaan, Gokhale & Zwaan 196S , 1972) 



is a well-known observational fact that implies the presence of a current sheet around the spot in the photospheric 



layers. Studies of sunspot structure 1989 indicate that a significant part of the currents is also concentrated into 



a current sheet in larger depths below the surface. The calculations presented above demonstrate how turbulent 
erosion of the surface of a magnetic flux tube can lead to the spontaneous formation of a current sheet around the 
tube. We thus propose that the sharp boundary of sunspots may be a natural consequence of turbulent erosion. 
As observations indicate that pores and sunspots tend to have well-defined boundaries from the beginning of their 
existence, this process of current sheet formation should occur below the surface, before the eruption of the tube. 

With the model developed in this paper we can attempt to understand some aspects of sunspot decay. That 
this process may be due to the eroding action of external motions was in fact already proposed by Simon & 



Leighton (1964). However, one must keep in mind that sunspots are very far from being slender flux tubes, so the 



application of the present models may only be suggestive, and the resulting formulae very approximative. 

The constancy of w implies a constant flux loss rate/tube surface area; therefore, if one introduces an artificial 
sink term in the diffusion equation outside the current sheet, thereby increasing the gradient and the flux loss, 
the decay rate should increase. This was indeed found in a numerical solution. A physical mechanism capable 
of such flux removal from the neighbourhood of a flux tube could be a transverse flow sweeping away the diffuse 
field (as with toroidal flux emerging through the convective zone). Again, such processes could only increase the 
effectiveness of turbulent erosion, thus the lifetime estimates given in these applications should be regarded as 
rather conservative upper estimates. On the other hand, such a removal of the outer diffuse flux may be regarded 
as a "renovation" of the initial state, so that the end of the period of effective renovation (i.e. the full emergence 
of the loop with verticalization of the legs) may be regarded as the instant of time corresponding to t = in our 
calculations. This may be approximately also the time when the spot area attains its maximal value, so "ro" may 
represent the maximal radius of the spot. 



Using equation (16) the lifetime of a spot in the limit of strong inhibition of turbulence is 



/w = 2 1 / 3 



Substituting here B ~ 3000 G, B e 



400 G, vo ~ 1000 km 2 /s (the granular value), one finds 
r /w = (r o /10 4 km) 2 x 10 d 



(31) 



(32) 



From this formula it is apparent that the total lifetime is proportional to the initial area of the spot, thus returning 
the well known linear area- lifetime relation ( pnevishcv 1938| ). The orders of magnitude are also consistent with 
observations: in particular, the total lifetime of large, recurrent spots of ~ 5 ■ 10 4 km diameter is found to be about 
2-3 months. 

The area-time decay curves resulting from a current sheet moving inwards with a constant velocity are neces- 
sarily parabolic, as can be seen in Figure ^ as well. This would run counter the common view that the sunspot 
decay proceeds linearly in time. There have been sugg estions of the parabolic nature of the sunspot decay laws on 
the basis of a statistical study of large data samples ( Moreno-Insertis fc Vazquez 1988 ). The scatter intrinsic to 
this kind of data, however, makes it difficult to distinguish a parabola with the curvature predicted above from a 
linear decay law (Martinez Pillct, Moreno-Insertis, fc Vazquez 1993 ). Nevertheless, in a recent analysis, Petrovay 
fc van Driel-Gesztelyi ( 1997 ) found decisive observational ev i dence in favour of a parabolic decay law. 

We note that in a parallel work Riidiger & Kitchatinov (1996) solved the 2- dimensional nonlinear diffusion 



equation. As far as it can be judged from their figures, their solution falls in the diffusive regime, with a nearly linear 
decay law (in apparent contradiction to the observations, cf. Petrovay & van Driel-Gesztelyi 1997). As however 
only results from a single run with one particular choice of initial conditions are available, and the uppermost 1500 
km of the convective zone (where B/B c would be highest) is not included in the model volume, it is presently not 
possible to judge the general validity of those findings. 

An interesting property of the erosion models is that the remaining lifetime of a sunspot with instantaneous 
radius res, 



or, equivalently, the area decay rate 



r cs /w = (r C s/ro)(ro/10 4 km) 2 x 12 d 



r CS (10 4 km) 5 

2-irrcsw = 7T — 

r 6 a 



(33) 



(34) 



does not only depend on res but also on ro, i.e. on the initial radius of the spot. This is caused by the fact that 
the value of w is fixed at an early stage of the decay, when the propagating solution is first adopted, and it does 



not change much thereafter. This implies that of two observed decaying sunspots with identical radii, the younger 
one (i.e. the one with a smaller maximal radius) should show a faster decay rate. One may speculate that this 
dependence on the original size may possibly be a contributing factor to the large intrinsic scatter in the decay 
rates of sunspots of a given size. 

From observations it is known that sunspots are surrounded by a radial outflow, called the moat. It is then 
interesting to investigate the effects of such an outflow on the solutions. The actual structure of the moat flow is 
not totally clear: since the field lines are inclined, the flow velocity could have a non-negligible component parallel 
to the field lines. However, only the component normal to the field may have a direct effect on its decay. We then 
use a more general form of the diffusion equation including an advection term due to the radial outflow of matter, 
and, within the one-dimensionality of the model, disregard any effects due to any other component of the flow. 
The resulting equation is: 



r "(B) f- 
or 



rB 



(35) 



which substitutes equation (Q). Observations sugg est that ou tside the penumbra the outflow has a velocity of 
~ 0.5-1 km/s, nearly independently of the radius. 197S , 1988 . On the other hand, the radial flow must have a 
sharp cutoff at the boundary of the strong field region. The function v r may then be represented as 



2v 



1 + (B/B c )<*» 

where vq is the outflow velocity at r e . Equation ( |l6| ) can be readily generalized to yield 
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(36) 



(37) 



We have obtained numerical solutions for this problem: the resulting field profile can be seen in Figure which 
presents the case for v — 1, a v — 5, a v — 3. These solutions confirm the accuracy of the estimate (0). Now the 
dimensionless value of v will be 0.5km/sxr /fo, i.e. crudely in the range 2 < VoTq/vq < 15. Consequently, a moat 
flow perpendicular to the field lines would imply an order-of-magnitude increase in the decay rate, in conflict with 
the observations. This leads us to conclude that, to the extent that the present models can teach us something 
about the process of sunspot de cay, the moat flow should be mostly parallel to the field lines. This conclusion also 
agrees with the observations of Skumanich, Lites, & Martinez Pillet (1994) that the apparent flux transport rate 
of the moat (in the form of moving magnetic features, MMFs) greatly exceeds the actual flux loss from the spot. 




Fig. 11. — Same as Fig. lc but including a radial outflow as in eq. (36), with vq = 1, a v = 5, Aio = 0.5. 



5. Conclusion 

An investigation by numerical and analytical means of the solution of the nonlinear axisymmetric diffusion 
equation ([!]) with v(B) functions that satisfy the minimum set of physical requirements listed in § |2.l| yielded the 
following main results. For sufficiently strong nonlinearity (i.e., for sufficiently strong reduction of v inside the 
tube) a current sheet is spontaneously formed around the tube within one diffusion timescale r\ / vq (ro is the initial 



radius of the tube, and vq is the kinematical value of the diffusivity) . The field profile in and inside the current 
sheet is well approximated by the analytical form-invariant propagating solution of the plane-parallel equivalent 
of equation ([j]) . This sheet propagates inwards with a velocity 

w ^ 2 -^^ (38) 

where Bq is the field strength just inside the current sheet, B e is the equipartition field strength. Accordingly, the 
lifetime of a tube with constant internal flux density is increased by a factor B / B e , independently of the value of 
the diffusivity inside the tube. 

On the basis of these results we performed approximate calculations of the magnetic flux loss from toroidal flux 
tubes lying at the bottom of the solar convective zone, and rising through the zone to the surface, respectively. 
It was found that the timescale of flux loss is ~ 1 month, comparable to the rise time, and very much shorter 
than the solar cycle. Consequently, toroidal flux bundles cannot be stored inside the convective zone proper for 
extended periods of time and the layer of flux storage must have turbulent diffusivities (and therefore turbulent 
velocities) significantly lower than the convectively unstable layer. Flux loops rising through the convective zone 
lose a significant fraction of their magnetic flux during the rise. 

While sunspots are far from being thin flux tubes, an application of our models to them still yields decay 
times comparable to those observed; besides, the linear area-lifetime relation is also returned. The inclusion of 
a moat flow with the observed velocities and perpendicular to the field lines however reduces the decay times to 
values incompatible with observations; thus, the moat flow may actually rather be more or less parallel to the field 
lines. Furthermore, turbulent erosion may explain the origin of the current sheet at the boundary of a sunspot. A 
curious feature of our model is that the decay rate does not only depend on the actual size of a spot but also on 
its original (i.e. maximal) size. This is the signature of a parabolic decay law only slightly departing from linear. 



Observations seem to confirm our prediction of a parabolic decay law (Petrovay & van Driel-Gesztelyi 1997). On 
the other hand, from the solutions obtained in this paper, we can estimate the Joule dissipation associated with 
the current sheet around a flux tube. Setting in the values corresponding to a sunspot, we obtain a photospheric 
power source that is 4 orders of magnitude lower than the solar radiative output on the corresponding area (an 
annulus of thickness 1" around the spot), far too weak to be observed. 

On the theoretical front, possible extensions of this work include a generalization by allowing the variation of 
physical quantities in the direction parallel to the tube, as well as tests of these results in three-dimensional MHD 
simulations. 

The authors are grateful to Dr Manfred Schussler for useful comments on this manuscript. This work was 
financed in part by the DGES project no. 95-0028 and by the OTKA under grant no. F012817. 
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